using System;
using System.Reflection;
using NBody;
using OptionParser;

//
// The code in this program is copied from John Dubinski's fakenfw program
//

[assembly: AssemblyTitle ("NBodyFakeNFW")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Calculate an NFW ellipsiodal system.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class FakeNFWOptions : OutputGetOptions
{
    [Option ("-n", "Number of particles")]
    public int N;
    
    [Option ("-m", "Cut-off mass")]
    public double Mcut;
    
    [Option ("-q1", "Axial ratio b/a")]
    public double Q1;
    
    [Option ("-q2", "Axial ratio c/a")]
    public double Q2;
    
    [Option ("-phi", "Longitude of ellipsoid minor axis")]
    public double Phi;
    
    [Option ("-inc", "Inclination of ellipsoid minor axis")]
    public double Inc;
    
    [Option ("-seed", "Random seed")]
    public int Seed;
    
    public FakeNFWOptions()
    {
        N = 100000;
        Mcut = 2.0;
        Q1 = 0.7;
        Q2 = 0.5;
        Phi = 60.0;
        Inc = 30.0;
        Seed = 0;
    }
}

public class NBodyFakeNFW
{
    public static void Main(string[] args)
    {
        FakeNFWOptions opts = new FakeNFWOptions();
        opts.ProcessArgs(args);
        
        // create a new system
        NBodySystem s = new NBodySystem(opts.N, 0.0);
        
        // cut radius
        double r1 = RadMass(opts.Mcut);
        Console.Error.WriteLine("Cut-off radius is {0}", r1);
        
        // set random seed
        Random random;
        if (opts.Seed == 0)
            random = new Random();
        else
            random = new Random(opts.Seed);
        
        double ci = Math.Cos(opts.Inc * Math.PI / 180.0);
        double si = Math.Sin(opts.Inc * Math.PI / 180.0);
        double cp = Math.Cos(opts.Phi * Math.PI / 180.0);
        double sp = Math.Sin(opts.Phi * Math.PI / 180.0);
        
        for (int i = 0; i < opts.N; i++)
        {
            double m = opts.Mcut * random.NextDouble();
            
            double rad = RadMass(m);
            
            double phi = 2.0 * Math.PI * random.NextDouble();
	
            double ct = 2.0 * (random.NextDouble() - 0.5);
            double st = Math.Sqrt(1.0 - ct * ct);
            
            double mass = opts.Mcut / (double)opts.N;
            
            double x0 = rad * st * Math.Cos(phi);
            double y0 = opts.Q1 * rad * st * Math.Sin(phi);
            double z0 = opts.Q2 * rad * ct;

            double x1 = ci * x0 + si * z0;
            double y1 = y0;
            double z1 = -si * x0 + ci * z0;
            
            double x2 = cp * x1 - sp * y1;
            double y2 = sp * x1 + cp * y1;
            double z2 = z1;
            
            Particle p = new Particle(mass, x2, y2, z2, x2, y2, z2);
            s.AddParticle(p);
        }
        
        s.AdjustForCentreOfMass();
        
        s.Write(opts.OutFile);
    }
    
    private static double RadMass(double m)
    {
        double eps = 1e-6;
		double dr = 1.0;
		double r1, r0 = 1.0;
		while (dr > eps)
        {
            r1 = r0 - (Math.Log(1.0 + r0) - r0 / (1.0 + r0) - m) / 
                      (r0 / (1.0 + r0) / (1.0 + r0));
            dr = Math.Abs(r1 - r0);
            r0 = r1;
        }
        
        return r0;
    }
}
